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PROBABILISTIC MICROMECHANICS FOR HIGH-TEMPERATURE COMPOSITES 
Grant Number: NAG-3-933; Period: Oct. 1, 1988-March 31, 1992 

Principal Investigator: J. N. Reddy 
Department of Engineering Science and Mechanics 
Virginia Polytechnic Institute and State University 
Blacksburg, Virginia 24061 

RESEARCH OBJECTIVE 

The three-year program of research had the following technical objectives: the 

development of probabilistic methods for micromechanics— based constitutive and failure 
models, application of the probabilistic methodology in the evaluation of various composite 
materials and simulation of expected uncertainties in unidirectional fiber composite 
properties, and influence of the uncertainties in composite properties on the structural 
response. The first year of research was devoted to the devdopment of probabilistic 
methodology for micromechanics models. The second year of research focused on the 
evaluation of the Chamis— Hopkins constitutive modd and Aboudi constitutive model using 
the methodology developed in the first year of research. The third year of research was 
devoted to the development of probabilistic finite dement analysis procedures for 
laminated composite plate and shell structures. 

TECHNICAL APPROACH AND RESEARCH ACCOMPLISHED 

To simulate variations and uncertainties in constituent properties and 
manu&cturing nonuniformities, as they would apply to high— temperature, metal— matrix 
composite materials and structures, the nonlinear behavior of the constituents was 
considered. The first phase of the research dealt with the linear micromechanics rdations 
for metal— matrix composites using reference temperatures and zero load values. 

The second phase involved the nonlinear characteristics (i.e. use of the powersaw 
rdations) and cyclic degradation effects. A Monte Carlo simulation routine has been built 
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around the METCAN program. The METCAN program (METal matrix Composite 
ANalyzer) uses the Chamis— Hopkins semi-empirical micromechanics equations in 
conjunction with a constituent based nonlinear power law to perform nonlinear analysis of 
fiber reinforced metal matrix composites. The research of the second phase was divided 
into three parts. The first part dealt with the linear studies, the second concerned with the 
nonlinear studies, and the third included comparisons between the Chamis— Hopkins 
micromechanics model and the Aboudi model. 

The third phase of research involved application of the research performed in the 
first two phases to the analysis of laminated composite plates and shells. A 3— D 
degenerated shell element formulation, incorporating the first-order shear deformable 
kinematics, is used in conjunction with the first-order second moment technique for 
probabilistic analysis. Linear elastic composite behavior is modeled as well as nonlinear 
material behavior using the micromechanics theory of Aboudi and the classical 
rate-4ndependent incremental plasticity, Random variables include: the ply stiffnesses, 
orientation angles, ply thicknesses, as well as the constituent micro properties such as the 
fiber and matrix stiffnesses and volume ratios. The probabilistic finite element model is 
used to quantify the variability and uncertainty in the response of composite shell 
structures. The computational procedure is then used to perform reliability analysis of 
composite shell structures. 

SIGNIFICANCE OF THE RESEARCH 

The probabilistic micromechanics approach and probabilistic finite element analysis 
computational procedures developed during this research are the first ones for laminated 
composite plates and shells. The procedures can be used in the analysis and design of high 
temperature, metal— matrix and ceramic composite materials and structures. The results of 
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this research play a crucial role in the establishment of increased performance and 
durability limits of high-temperature composite structures. 
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TECHNICAL DISCUSSION 

1. Intiodnction 

The current emphasis for computational analysis of high temperature metal matrix 
composites involves micromechanics theories with nonlinear behavior occurring at the 
constituent level. Interphase regions are included in these models, as studies have shown 
that residual stresses that occur in the composites can be controlled by appropriate choice 
of an interphase material and thickness. Both ply and constituent level failure theories are 
employed in the resultant schemes. The procedures, due to the use of constituent level 
material properties as the primitive variables, can be used to predict the qualitative 
behavior of new composites prior to any experimental measurements. Material properties, 
strengths, and stress-strain data are useful predictions. However, since uncertainties are 
present in the constituent properties, volume ratios, and fabrication variables, the results 
lack the variability that is always present in experimental data (see [1—4]). 

A previous study by Stock [4] proposed that constituent level uncertainties be 
modeled using statistical distributions and the resultant composite micromechanics modds 
and laminate theories be used to predict the statistical scatter in experimentally observed 
composite properties. A computer code developed at NASA Lewis, called ICAN 
(Integrated Composite Analyzer), for comprehensive linear analysis of multilevd fiber 
composite structures was the deterministic basis of the simulation modd. Monte Carlo 
simulation methods were used to model the probabilistic distributions in the primitive 
variables, and distribution results for material properties and strengths of a graphite— epoxy 
ply were presented. The probabilistic micromechanics model proved quite useful in 
predicting the experimental scatter of these linear composite structures. 

Advanced composites such as the metal matrix composites are used in aerospace 
structures where operating temperatures exceed the limits of organic matrix composites. 
Therefore, researchers have focused on the devdopment of suitable constitutive theories for 
these materials. Dvorak and Bahei-El-Din [5-9] introduced a computationally simple 
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niicToni6cli3iiic8 iiiodsl b&sed on. what they refer to as a ’'vanishing fiber diameter model”. 
The model has been used in the analysis of a variety of laminated composites [10—12]. 

The power-law model of Hopkins and Chamis [1] is based on the assumption of 
strain compatibility and stress equilibrium. In this model the representative volume 
element is subdivided into several distinct phases. The different phases are used to 
represent the fiber, an interphase region and three matrix subregions. The matrix 
subregions were introduced to characterize the through— the— thickness variation of the 
matrix phase. 

The "periodic hexagonal array" (PHA) model of Teply and Dvorak [13,14] is 
developed by considering a hexagonal fiber array as the representative ceU. This leads to 
the use of a triangular representative volume element which is analyzed utilizing the finite 
element method. Both displacement and equilibrium formulations are examined to 
establish upper and lower bounds on the instantaneous properties of the composite. 

Aboudi [15-20] developed a comprehensive micromechanics theory which has its 
origins in the work of Achenbach and his colleagues [21,22]. The theory is applicable to 
several t 3 q>es of composites, such as particulate, short— fiber or continuous— fiber 
reinforcement. Unlike most micromechanics theories, the Aboudi theory is not restricted 
to the assumption of perfect fiber/matrix bonding. Instead, an infinitely thin elastic 
interface is assumed which allows the simulation of fiber/matrix damage or an interphase 
region. The Aboudi micromechanics theory can be viewed as a variational formulation of 
the displacement or stress boundary value problem at the micromechanics level (see 
[23-25]). 

Another approach taken by many researchers when modeling metal matrix 
composites does not involve micromechanics level analysis. They attempt to model the 
material nonlinearity using macroscopic yield criteria to account for the amsotropic 
behavior. Many authors have made contributions in this area; however, only a few are 
discussed here for brevity. Hill’s anisotropic theory of plasticity [26] has received much 
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attention. This theory was based on a generation of the von mises yield criterion which 
assumed yielding was independent of hydrostatic stress and that plastic flow was 
incompressible. Whereas these assumptions are standard and experimentally validated for 
metals, they have been shown to be incorrect for some composite materials. A recent yield 
function was introduced by Sun et al. [27,28], which was selected for this work. The 
function is quadratic in stresses and, in general, excludes the assumptions of 
incompressibility of plastic strains and that no yielding is caused by hydrostatic stresses. 

It is the objective of this study to quantify the variability present in metal matrix 
composites. In the first model, the probabilistic micromechanics concept developed by 
Stock [4] to apply to the nonlinear constitutive behavior of high temperatme metal matrix 
composites. The Monte Carlo procedure is employed in conjunction with the 
METCAN program (METal matrix Composite ANalyzer) developed at NASA Lewis 
Research Center. Assumed statistical distributions for material properties, volume ratios, 
and nonlinear parameters thus allow prediction of experimental scatter in metal matrix 
composite behavior. In the second model, a macromechanics (orthotropic) elastoplastidty 
theory is combined with a continuum shell element (see Reddy and his colleagues [29,30]) 
and a probabilistic methodology is incorporated to account for uncertainties in material 
properties. The resulting formulation is used to quantify the variability of the structural 
response of materially nonlinear composite laminates. Numerical examples based on the 
formulation are also discussed here. 

2. Miciomechanics Constitutive Model 
2.1 Micromechanics Formulation 

The deterministic part of the model consists of the computational scheme contained 
in the METCAN program. It employs the Chamis— Hopkins semi— empirical 
micromechanics equations [1] in conjunction with a constituent based multi— factor 
nonlinear material model [2,3] to perform nonlinear analysis of fiber reinforced metal 
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matrix composites. As shown in Fig. 1, a closed-^oop analysis logic is used in which 
nonlinearities are applied at the micromechamcs (constituent) level, ply properties are 
constructed in the synthesis stage, laminate theory is used to build laimnate properties, 
global analysis is performed, and the laminate is decomposed again to the constituent level. 
Note that the Monte Carlo sampling is performed at the constituent level. In the 
decomposition stage, the ply stresses are computed and used to calculate the stresses in the 
constituents (microstresses). The constituent material properties are updated based on 
their dependency on various factors in the nonlinear material model including stress, 
temperature, time, mechanical and thermal cycles. The resultant nonlinear analysis 
involves incremental loading and time stepping in which iterative convergence is checked at 
the micro, ply and laminate scales. 

The multi-factor material model is presented in Figure 2 along with the 
representative subcell of the micromechanics relations. This material model has been 
selected to consistently represent the in— situ behavior of all constituent properties. More 
details on this model can be found in the work by Hopkins [2]. The unit cell presented in 
Figure 2 illustrates the various subregions assumed in the micromechanics theory to 
capture the effects of matrix, interphase, and fiber interaction. 

2.2 Monte Carlo Simulation 

Monte Carlo simulation is a computational technique often used to simulate random 
processes. In a general sense, it is defined as any computer simulation involving random 
numbers for solving stochastic problems. In the procedure, the computer is used to 
generate independent statistical samples for each random variable, which are then fed into 
the model. Each sample can be thought of as an independent deterministic experiment, 
which is processed by the model to yield the results of the experiment. Each sample is 
drawn from an appropriate probability distribution. Following the simulation process, the 
output data is statistically analyzed to estimate the true characteristics of the model. The 



10 


law of large numbers assures that if sufficiently large sample sizes are taken, the results will 
converge to the true population statistics of the model. 

3. ’M'arTQTTierlia.iii Qt Constitutive/Structural Modd 
S.l Orthotropic Plasticity Formulation 

In order to model orthotropic dastoplastic behavior from a macroscopic viewpoint, 
a quadratic yield function introduced by Sun, et al. [27,28] is adopted. The yield function 
is quadratic in the stresses and employs the associative flow rule and isotropic hardening. 
The plane stress radial return algorithm of Simo and his colleagues [31,32] was modified to 
include this new yield function. This algorithm was chosen due to its accuracy, improved 
global convergence rates, and compatibibty with the probabilistic finite element routines. 


Governing Equations 

According to Sun [27], the orthotropic yield function is given by 


^^ 11^11 + h 2^22 ^ 33^33 ^^ 2 ^ 11*^22 ^^ 3 *^ 11^33 


2 2 2 " 
+ ^^3 ^ 22^^33 ^^ 44^^23 ^^ 55^^13 ^^ 66 ^^ 12 ^ “ ^ 


( 1 ) 


where Y is a state variable and crij are the stresses in principal material coordinates. The 
coefficients aij are constants, which are determined from experimental data and control the 
amount of anisotropy in the plastic behavior. This yield criterion does not include the 
assumption of incompressibility of plastic strains or that hydrostatic stresses result in no 
yielding or plastic deformation. The function also reduces to the von Mises criterion or the 
Hill yidd function for orthotropic materials [26] with appropriate sdection of the aij 


values. 
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T1i 6 as 80 ciativ 6 flow rule assumption allows the incremental plastic strains to be 
stated as 




ij 


( 2 ) 


The shell finite element employed is a continuum element developed from the 3— D 
elasticity theory using the kinematics of the first order shear deformation theory (see 
[29,30]). In other words, transverse normals remain straight and inextensible, so that £33 = 
0. Consequently, <733 does not enter the formulation (because the strain energy of the shell 
due to £33 is zero). 

In vector notation, the stress and strain tensors can be written as 


{< 7 ^} - {^11 <^22 *’^12 *^13 ‘’' 23 ^ 

{ c } = {^11 ^22 ^^2 ^^13 ^^23^^ ’ ^ ^^11 ^22 ^^12 ^^13 ^^23^ 

The components of the back stress qy (included to model kinematic hardening), and the 
relative stress rm = trij - qij are expressed in the vector form as 

{q} = {Qh fi22 ^112 *^13 ^^ 23 ^ > {^} = ^22 ^^12 ^13 

Thus the governing elastoplastic equations in vector form can be expressed as (see [31]) 

{£} = {£*} + {C") 

{<r} = lQ]{c*} 

{i'*} = MPHO} 

{q} = 7|H'[Pj{q} 
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£ = ^’iFip|W-jyV)<o 
“=7l|{»f[PlW!‘/^ (5) 

where 7 is the time derivative of the plastic load parameter, H' is the kinematic hardening 
modulus, and the parameter Y represents the hardening law in terms of the equivalent 
plastic strain a. The matrix [Q] is the elastic constitutive matrix, adjusted for the 
constraint <733 = 0. The effective stress a can be expressed as 



and matrix [P] is given by 

^11 0 0 0 

a-t n ^00 0 0 0 

0 0 2agg 0 0 (7) 

000 2ag5 0 

0 0 0 0 2a^4 . 

Loading and unloading conditions are stated in the Kuhn— Tucker form [32] by 
requiring that 

f<0,7>0,^ = 0 (8) 

For an elastic process, f < 0 and 7 = 0. For a plastic process, we have f = 0 and 7 > 0. 
These two conditions are generally valid, for a loading or unloading state. 

Incremental Formulation 

The ordinary differential equations of time in Eq. (5) can be numerically integrated 
using a backward Euler difference scheme over the time interval (tn,tn+i). Letting 7n+i = 

7n*iAt (the plastic load parameter) and I = J the strain at tn*i can be 
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written in terms of the strain at tn and the gradient of the inaemental displacements, 

where "V" denotes the differential operator used in the definition of strains. A trial stress 
state is assumed by freezing plastic flow so that the entire step is purely elastic. The trial 
stress then becomes 

{'Ll> = + \+lPK’>n+l> 

"''n+l 

“n+1 “n + '^B+l^n+r 

In order to perform the incremental updates required in equations (10), the plastic 
load parameter (Lagrange multiplier) 7 must be determined. It is found by enforcing the 
consistency condition at time tn*i, i-e., 


= I ^n+1 ^ 


( 11 ) 


The actual hardening functions used in the program are those recommended by Simo and 
Hughes [31], 

Y(a) = iMa + ay + (K^ - “ exp(-Aa)] (12) 


and by Sun [27], 
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Y(a) = s[a+ 

where H is the hardening modulus, the uniaxial yield stress, and A are other input 
parameters, and 

H(a) = (l-)9)Ha. 

Here 0 denotes the fraction of kinematic and isotropic hardening desired, i.e. /3 = 1 denotes 
purely isotropic hardening, and ^ = 0 denotes purely kinematic hardening. Equation (11) 
is solved at each gauss point in the structure for 7 by a local Newton iteration procedure, 
as it is generally a nonlinear scalar equation. 

The global equilibrium is obtained by using Newton— Raphson iteration. This 
requires that the tangent moduli be known in the form 

Simo Hughes [31] developed tangent moduli, which are consistent with the integration 
algorithm previously discussed. For finite values of load step size, they showed that the 
consistent elastoplastic tangent moduli preserved the quadratic rate of asymptotic 
convergence that is characteristic of Newton’s method. Use of the continuum tangent 
moduli derived independent of the algorithm loses this convergence rate. 

S.2 Continuum-Based Finite Element 

The incremental equations of a continuous medium are formulated based on the 
principle of virtual displacements and the total Lagrangian description. The detailed 
description can be found in [29,30]. The final incremental equilibrium equations for an 
element are given by 
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aKL] + [KNLlH^> = W-W 

where {A} is the vector of nodal incremental displacements in an element, and [K^^], 
and {F} are defined by, 


IKJ = LlB/WlPlJdV 

j y 

{F} = [ [BJ{S}dV. (15b) 

In the above equations, [B^^] and are linear and nonlinear strain-displacement 

transformation matrices, [Q] is the incremental stress-strain material property matrix, {S} 
is a vector of 2nd-Piola-Kirchho£f stresses, and {R} is the external load vector. All 
matrix elements refer to the deformed state with respect to the original undeformed 
configuration. 

After assembly, the following linearized versions of the actual equilibrium equations 
are obtained 

[K]{A} = {R}. (16) 

where 

[K] = [Kj,] + [KniI 

[R] = {R} - (F) (17) 

The Newton— Raphson method is employed to solve the linearized equations iteratively 
iint.i1 the actual equations of motion are satisfied to a required tolerance. 
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In the process of evaluating the integrals in equation (15), Gauss quadrature is used 
in the membrane directions of the shell, and in the thickness direction when plasticity 
occurs. When material behavior is elastic at a membrane location through the entire 
laminate thickness, then explicit integration is performed in the thickness direction. 

S.S Second-Moment Formulation 

The second— moment perturbation method as developed by Liu, et al. [33,34] is 
summarized in this section for geometric and material nonlinear time independent 
behavior. 

The equilibrium equation for the probabilistic finite element model can be expressed 
as 

{F({A},{b})} = {R({b})} (18) 

where {F} is the internal force vector, {A} is the displacement vector, {R} is the external 
force, and {b} is a discretized vector of the random function b(x), where x is a spatial 

coordinate {x}. The random function b(x) is expanded using shape functions ^i(x): 

b(x)= S V'jWbj (19) 

i=l * ■ ‘ 

where bi are the nodal values of b(x). Generally the random quantity b can be a material 
property, geometric dimension, or a load. 

The probabilistic finite element analysis is carried out by applying the 
second— moment method. The vectors in equation (18) are expanded about the mean of the 
random function b via Taylor’s series, and the following perturbation equations are 
obtained: 
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Zeroth-order equation: 

{F} = {R} 

(20) 

First-order equations: 

“ {^}b. ~ ~ l,...,n 

(21) 


where the bar over the symbols represent evaluation of that term at the mean value of b, 
qnH the subscript bi represents the derivative with respect to bi- Also [K ] is defined to be 
the tangent sti&ess matrix, 


[kT]=||^ (22) 

Once {A} and {A}b. are obtained by solving equations (20) and (21), the mean and 
autocovariance matrices for the nodal displacements can be determined (see [38]). 

The solution procedure involves a consecutive solution of equations (20) and (21). 
After the deterministic zeroth— order equation (21) is solved, either once for the linear case 
or iteratively at each load step in the nonlinear case, the generalized displacement vector 
{A} is used to perform the perturbation solutions in equation (21). In (21), there are as 
many solutions for {A}b. as there are number of nodes in the model. In addition, the 
computations must be performed for each layer in the model, as a particular random 
function is assumed to be independent from layer to layer. If there are n nodes and P 
layers in the model, (n + 1)P more matrix solutions are required at each load step for each 
random variable. This is not as expensive as it seems, because the stiffiiess matrix [K ] is 
inverted once and used in (21). 

For materially nonlinear problems the residual vector can be functionally 
represented in the form 

{F} = {F({A},{b})} (23) 

Temporarily dropping the vector braces, and difierentiating F with respect to b, we obtain 
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dF _dF ,dF dA 

3F W 


(24) 


The first term on the right hand side involves the explidt derivative of F with respect to b. 
This is the derivative that must be evaluated on the right side of equation (21). This 
derivative can be expressed exactly when elastic behavior exists. However, when material 
nonlinear behavior is present, finite difference derivatives are used (as in [33]) once again 
due to the complex nonlinear relationship between stress and strain. Since orthotropic 
plastidty is present in certain layers, the radial return algorithm discussed previously is 
considered to be very good for this purpose (see [33]). This is true since at a particular 
load step, the plastic strain and effective plastic strain are "frozen” and only updated after 
equilibrium conversion is achieved. This update can be done after the finite difference 
calculations are made, so that the true effect of perturbing the random variable about its 
mean is measured. Other advantages of the radial return method are the increased 
accuracy involved in the stress recovery routine and the algorithmic compatible tangent 
moduli which results in a more accurate tangent stiffiiess matrix. This accuracy is 
important in computing the perturbation derivatives [e.g., {A}b. in equation (21)]. 

Sources of randonmess in the model can be material properties, geometric 
dimensions, or loads. In the present study, the ply thickness, ply angle and material 
properties are selected as the random variables. The random material variables include the 
engineering properties Eii,E 22 j^^ 12 >Gi 2 |Gi 3 >G 23 > nnd the plasticity parameters (Xy and H 
(uniaxial yield stress and hardening modulus). The loading is considered to be 
deterministic throughout this study. For details regarding the incorporation of these 
variables, see Reference 38. 

The probabilistic finite element procedure developed herein has the ability to model 
the spatial correlation involved in random fields. This technique is discussed in detail in 
[38], where examples involving random fields were included. However in the present work 
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it is assumed that fully correlated fields exist in each ply which degenerate to single 
random variables, and each random variable is independent from ply to ply. 

4. Applications 

4.1 Micromechanics Constitutive Model 
Linear Analysis 

First, simulations were performed in which the nonlinear multi— factor material 
model was not activated. As in Stock’s work [4] histograms, cumulative distribution, and 
(vtnfii^pTirp interval curves were used as methods of illustrating various characteristics of a 
single ply. PlOO high modulus graphite fiber and a copper matrix were selected as the 
constituents for this linear work. In general, sample sizes of 50 were chosen as statistically 
significant to show the trends of interest and remain economical. 

Histograms and cumulative distribution plots were made for the two cases given in 
Table 1. Case 1 represents a narrow range distribution of the constituent properties, and 
Case 2 a wide one. A deterministic case was run using the mean values from Cases 1 and 2 
(given by either y.ot Sample results are given in Figures 3a-d for the in— plane shear 
modulus G ^2 results. These are compared to the deterministic case result of 32.2 GPa to 
assess the sensitivity of the ply properties to changes in the distribution parameters. 

Confidence interval curves are used to investigate and illustrate various effects on 
the resultant ply properties. These effects include fiber and matrix strengths, fiber 
orientation, fiber stiffnesses, and interphase thickness. Some examples of these effects are 
presented next. 

Fiber Strength Effect . Using the properties of Case 1 as the base, the shape 
parameter of the WeibuU distribution for fiber strength was perturbed to show the effects 
on the ply properties. The fiber volume ratio was kept deterministic and the simulation 
was performed for a range of these values. Figure 4 contains representative results for the 
ply longitudinal tensile strength. The solid lines represent the means of the samples, and 
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the stars above and below represent the upper and lower bounds of the 95% confidence 
interval estimates. 

Matrix Strength Effect . Similar variations are made on the matrix strength shape 
parameters with the intent of showing the effects of these changes on ply strengths. These 
results are given in Figures 5 through 7. 

Fiber Orientation Effect . Three values of the distribution parameter on fiber angle 
were simulated to determine the effect on the ply properties. These curves are given in 
Figures 8 through 11. 

In general, from reviewing the previous confidence interval plots, the metal matrix 
composite ply properties are not very sensitive to fiber misalignment due to the closeness of 
the matrix modulus to the fiber modulus. It is also easy to see that ply longitudinal tensile 
strength is significantly affected by the shape parameter of fiber strength. Without further 
detail, it is evident that much information concerning the variability of this metal matrix 
composite has been quantified. The next step is to determine effects of allowing the 
nonlinear material model to become active. 

Nonlinear Analysis 

For any given particular ply property, the governing equation is the micromechanics 
relation shown in Fig. 2. This equation relates the ply property to the constituent 
properties, i.e., the fiber, matrix, and interphase properties. These constituent properties 
are in turn governed by the nonlinear multi— factor material model (power law), thereby 
making them a function of the local temperatures, stresses, stress rates, and mechanical 
and thermal cycles (see Figure 2). The constituent power laws contain certain allowable 
parameters such as melting temperature, strengths, allowable stress rates and cyclic 
thermal and mechanical strengths. In addition, the individual terms in the power law are 
raised to an exponent or power. 
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In this section, in addition to the random variables used in the linear analyses, the 
allowable parameters in the power law terms and the exponents were included as random 
variables. The material studied is SiUcon Carbide fiber (SCS- 6 ) in a Titanium Aluminide 
matrix (TI15). Table 2 contains information on the distribution parameters for each 
random variable. For illustrative purposes, confidence interval plots are used again to 
investigate the effects on the ply longitudinal tensile strength. The fiber volume ratio is 
deterministic and a sample size of 50 is again chosen. 

Longitudinal Tpnsilft Strength . The micromechanics equation for longitudinal 
tensile strength is given below for ease of discussion: 
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The symbols K, 7 and D indicate volume ratio, an empirical constant, and fiber diameter 
respectivdy. The nonlinear power law rdationship for the appropriate constituent 
properties are. 
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Similar relations hold for Edn and Efn. In order to have a basis for comparison. Figure 12, 
which contains no nonlinear effects, is included. 

First the effects of the constituent property Sfut are investigated. In Figures 13 and 
14, the nonlinear temperature term is activated by elevating the ply to 811* K temperature. 
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Here the coefficient of variation (COV) of the temperature exponent n was varied from 1% 
to 30%, and the COV of the fiber melting temperature was varied from 1% to 10%. One 
can see that generally the values were reduced when compared to the linear curve; 
however, little difference occurred for the perturbation on the melting temperature or the 
exponent n. As for the stress nonlinear power term, the Weibull distribution parameter for 
Sfiit was varied between 20 and 10 (narrow and wide) when a ply uniaxial stress of 68.9 
MPa existed (see Figure 15). In order for the fiber tensile strength to behave linearly to 
stress changes as expected, the exponent m in the micromechanics equation was chosen to 
be zero (see Fig. 2). Thus the differences in the two curves shown in Figure 15 are simply 
due to the influence of the micromechanics equation and not the power law equation. 

Next, the nonlinear behavior of the constituent Enii was studied. Figure 16 shows 
the effect of simply varying the COV of Enu without the power law terms being activated. 
One can see that the effects of changing Emu’s distribution parameters are small. Figures 
17 and 18 contain the results of activating the temperature term and perturbing the 
exponent n and the matrix melting temperature COV’s in equation (27). Little changes 
due to the perturbed quantities are noticed for these results. 

Figure 19 is an example of activating the mechanical cycle nonlinear power term. In 
this plot the Weibull distribution parameter for the mechanical cyclic strength was 
varied from 20 to 10. In this example, 5 x 10® cycles (50% of the assumed strength) was 
used. From this result, it can be concluded that the probabilistic distribution of cyclic 
strength is significant for this high degree of nonlinearity, as would be expected. 

Obviously, the possibilities for analytical examples and investigations are 
voluminous. The previous examples are meant only to be representative of the capability. 
In addition, as in Stock [4], the individual ply could be broken down into subplies and the 
simulation procedure performed at this scale. Laminate theory can be used to build the ply 
properties, and in this way uncertainties caused by many fibers per layer or variations 
through the thickness is effectively modeled. 
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J^.Z Macromechanics Constitutive! Structural Model 
ARALL Laminate Tension Specimen with Hole 

ARALL laminates are high strength hybrid composites devdoped by Alcoa for 
aerospace applications. Figure 20 illustrates the architecture of the laminate, which 
consists of thin sheets of high strength aluminum alloys bonded to high strength aramid 
fibers using a special epoxy resin. The ARALL launinates have increased fatigue life and 
fatigue crack growth properties over monolithic aluminum; the outer aluminum layers 
provide impact Ha.Tn a.gft and moisture protection that would be a problem for typical fiber 
composite materials. In addition, increases in strength and lower doisities are achieved as 
compared to monolithic aluminum [38]. Here we analyze the ARALL laminates using the 
macromechanics elastoplastic formulation discussed earlier. 

Figure 20 also shows the description of the analytical modd used to represent the 
stifihess of the ARALL laminates. The aramid epoxy layers are divided into fiber— rich and 
resin— rich layers. Table 3 contains the properties and statistical distributions for the 
aluminum, aramid epoxy fiber— rich, and aramid epoxy resin— rich layers. Experimental 
tension test results [37] are compared with the analytical results in Figure 21. From the 
figure it is observed that the aramid epoxy behavior is linear and the analytical linear 
comparison is very good. The 7075— T6(L) aluminum behavior is dastic perfectly— plastic 
and the analytical modd Avith a yidd stress of 78 ksi agrees very well except near the point 
of first yidd. As for the ARALL-1 results (-1 indicates 7075-T6(L) aluminum is used) 
the analytical model with ideal plastidty for the aluminum layers and linear dastic 
aramid epoxy layers generally exhibits the same behavior as the experimental results 
except the 0 degree laminate analytical modd underpredicts the stifihess after yidd and 
the 90 degree laminate modd overpredicts the stiffness after yidd. For the purpose of this 
example the analytical model is considered acceptable and will be used to study the mean 
and variance response of an ARALL tension specimen vdth a hole. 
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Figure 22a shows the finite element model and dimensions of the tension specimen 
with a hole problem. The same material properties and material model from the previous 
discussion are used in this problem. The probabilistic analysis assumed a fully correlated 
random field for each random function in each layer. Figure 23 contains the mean and 
standard deviation of the longitudinal Cyy strain at the hole edge (point A) for the case 
where all aramid epoxy layers are aligned at 90 d^rees to the load. The figure also shows 
the breakdown of standard deviations for all the significant random variables. Since the 
fibers are oriented at 90 degrees to the load and to the strain Cyy, then the aluminum 
properties tend to dominate. It is interesting to note that even though no bending occurs 
in this problem, the ply thickness of the aluminum layers is dominant after yield. The 
alunoinum yield stress and elastic modulus are also important. Figure 24 contains similar 
results except now the fibers are aligned with the loading direction. While the aluminum 
yield stress and ply thickness random variables are still significant, the aramid En and ply 
t hicViK^s random variables are now equally important. These results illustrate the role the 
individual random variables play in the total variability of this type of ARALL structure. 

Boronf Aluminum Tension Specimen with Hole 

A Boron/ Aluminum laminate was selected to illustrate the use of the macroscopic 
orthotropic plasticity formulation. The same problem dimensions (except for thickness) 
were used as in the last example, however, a different mesh was used that placed gauss 
points along the x— axis (see Fig. 22b). Rizzi, et al. [28] conducted an experimental and 
analytical study of this specimen, and provided experimental measurements for the 
orthotropic elastic constants as well as the ay values in the yield criterion and the 
haTHpning parameters in the isotropic work hardening model. These values are all stated in 
Table 4 and are used in the present analytical model. It should be noted that the ay 
values used in this study differ from those given in the reference by a factor of 2/3 due to a 
minor difference in the formulations. Figure 25 contains a comparison of the analytical 
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results from the present study and experimental results from [28] for the longitudinal strain 
cyy along a radial line (x-axis) 90 degrees to the loading. The agreement is slightly worse 
than that obtained in [28], but is probably due to the difference in element formulations 
and the classical incremental plastic stress routine used versus the radial return algorithm 
used here. Yielding occurs after 1000 lbs, and the agreement worsens as the loading is 
increased. However, the results are still considered quite good. 

Using the random variable statistics stated in Table 4, the first— order 
second— moment probabilistic method was used to evaluate the mean and variance of the 
Cyy strain response. Once again the probabilistic analysis assumed a fully correlated 
random field for each random function. Figure 26 shows the analytical mean €yy strain for 
the 2500 lb and 1000 lb load values with the plus or minus one standard deviation points 
included. It is obvious that the sensitivity of Cyy to the random variables increases both 
with the load and as the location moves closer to the hole. Figure 27 contains a plot of 
both the mean and standard deviation of the Cyy strain at the location A on the model 
versus load. The breakdown for each random variable is presented as well. Since only a 
single layer is used, then the ply thickness could not be considered a variable here. The 
most significant random variable is the plastic hardening modulus H, with E 22 and the 
yield stress important as well. Note that E 22 is significant since the fibers are 90 degrees to 
both the loading direction and to Cyy. This example can be extended to include the aij 
plastic yield coefficients and the hardening parameter A as random variables since they are 
also experimentally measured quantities with uncertainties. 

5. Summary 

A probabilistic analysis procedure for constitutive behavior of metal matrix 
composites based on the METCAN program is developed. The procedure can be used to 
simulate manufacturing nonuniformities and uncertainties in constituent properties to 
quantify their overall effects on the composite. Studies involving both linear and nonlinear 



26 


effects on the thermoelastic and strength properties of two different metal matrix 
composites were performed. For the case of linear behavior the contributing constituent 
variations were constrained to the framework of the micromechanics model. Thus, cause 
and effects for the hnear behavior were easy to demonstrate, so that the relative 
importance of each material variable could be identified. As for the nonlinear effects, since 
the constituent-based nonlinear material model became active, variations were induced not 
only by the probabilistic distributions of the constituent properties but also by the 
distributions of the nonlinear power term parameters such as the melting temperatures and 
exponents. Thus the nonlinear behavior was really a blend of the variations in the 
micromechanics model variables and the nonlinear power law variables. It is easy to see 
how this procedure could be used to aid in material characterization and selection to 
precede and aid in expenmental studies. Much of the results presented have been based on 

assumed distributions, and thus are intended to be examples illustrating the power of the 
method. 

A formulation based on a macromechanics orthotropic elastoplasticity theory is also 
presented. A nonlinear probabilistic finite element analysis procedure including 
dastoplastic constitutive behavior is developed. The first-order second-moment method 
for probabilistic finite element analysis was combined with a continuum shell element 
which includes the effects of shear deformation. Both ARALL and Boron/ Aluminum 
plasticity problems were investigated, and the variability of these composites was 
quantified for a tension specimen with a hole. 
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TABLE 1 

Inpnt Statistical Parameters for Graphite Copper 


Iiipiit Case_l C3se_2 


Normally Distributed Variables 


cr 

a 

Ply Angle (degrees) 

0.0 

5.0 

10.0 

Fiber Volume Ratio (FVR) 

0.5 

0.1 

0.2 

Ejjj(GPa) 

723.9 

36.2 

72.4 

Eq2(GP») 

6.2 

0.3 

0.6 

Gfi2(GPa) 

7.6 

0.4 

0.8 

G^jCGPa) 

4.8 

0.24 

0.48 

Ej^(GPa) 

122.0 

6.1 

12.2 

E^(GPa) 

275.8 

13.8 

27.6 

Interphase % of fiber diameter 

0.10 

0.005 

0.01 

WeibuU Distributed Variables 

P 

A 

A 


2240.8 

20 

10 


1378.9 

20 

10 


172.4 

20 

10 


172.4 

20 

10 


172.4 

20 

10 


86.2 

20 

10 

SJMPa) 

220.6 

20 

10 

S^(MP®) 

131.0 

20 

10 

Sj(MPa) 

103.4 

20 

10 

Sd,(MP®) 

68.9 

20 

10 

Gamma Distributed Variables 

P 

A 

A 

Void Volume Ratio (VVR) 

0.33 

3.0 

5a0 
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TABLE 2 


Input Statistical Parameters for SCS— 6 TI15 


Normally Distributed Variables 
Ply Angle (degrees) 

%l(GPa) 

Ej2j(GP») 

G^(GPa) 

Em(GPa) 

Ej(GPa) 

TmiCK) 

TMdCK) 

Fiber Exponents 
Matrix Exponents 
Interphase Exponents 
Interphase % of fiber diameter 

“fll 

“{22 

“m 

“d 

Wpihnll Distributed Variables 
Sm,(MPa) 

S(ii^(MPa) 

S(j2s(MPa) 

S^(MPa) 

S^(MPa) 

Sj(MPa) 

S^,(MPa) 

rtamma Distributed Variables 


U 

0.0 

349.6 

349.6 

146.9 

146.9 

84.8 

275.8 

2755.4 

1255.4 

2199.9 
0.25 
0.50 
0.50 
0.10 

0.12E-5 

0.12E-^ 

0.45E-n5 

0.5E-5 

& 

3350.9 
3350.9 
3350.9 
3350.9 

1675.4 
1675.4 

896.3 

627.4 

103.4 
68.9 

1.0E6 

0.33 


0.5 

17.5 

17.5 

7.3 

7.3 

4.2 

13.8 

107.2 

32.2 

79.4 

0.0125 

0.025 

0.025 

0.005 

0.60E-7 

0.60E-7 

0.225E-6 

0.25E-6 

A 

20 

20 

20 

20 

20 

20 

20 

20 

20 

20 

20 

A 

3.0 


Void Volume Ratio (VVR) 



32 


TABLE 3 

Material Properties and Statistics for ARALL— 1 Laminate 

Constituents 


Random 

Variable 

Mean 

Standard 

Deyiation 

Coefficient of 
Variation 

Aluminum (7075— T6L) 



E 

10.4x10° 

5.2x10° 

0.05 

u 

0.3 

1.5x10"^ 

0.05 

* 

ay 

7.8x10^ 

3.9x10^ 

0.05 

* 

6 

1.2x10“^ 

6.0x10"^ 

0.05 

Aramid Etx)xv fiber— rich layers 



^11 

12.549x10° 

6.2745x10® 

0.05 

^22 

0.76525x10® 

3.82625x10'^ 

0.05 

^12 

0.28955x10® 

1.44775x10^ 

0.05 

*^12 

0.3458 

1.729x10”^ 

0.05 

^13 

0.28955x10® 

1.44775x10^ 

0.05 

♦ 




6 

0* ,90- 

2" 

— 

♦ 

S 

5.6x10“^ 

2.8x10“^ 

0.05 

Aramid Eooxv resin— rich layers 



Ell 

2.1972x10° 

1.0986x10® 

0.05 

®22 

0.48219x10® 

2.41095x10'^ 

0.05 

®12 

0.15717x10® 

7.8585x10® 

0.05 

*^12 

0.3749 

1.8745x10“® 

0.05 

°13 

0.15717x10® 

7.8585x10® 

0.05 

^23 

0.15576x10® 

7.7880x10® 

0.05 

* 

9 

0° ,90* 

2* 

— 

♦ 

6 

1.416x10”^ 

7.08x10“^ 

0.05 


*ay indicates yield stress, 0 indicates fiber orientation angle, and 6 

indicates ply thickness. 

Units are in psi and inches where appropriate. 


33 


TABLE 4 

Material Properties and Statistics for Boron/ Alnminnm 

Laminate 


Random 

Variable 

Mean 

Standard 

Deviation 

Coefficient of 
Variation 

Ell 

29.4x10® 

1.47x10® 

0.05 

®22 

19.1x10® 

9.55x10® 

0.05 

Gi2 

7.49x10® 

3.745x10^ 

0.05 

I/- ^ 

0.169 

8.45x10“® 

0.05 

12 




^13 

7.49x10® 

3.745x10 

0.05 

^23 

7.49x10® 

3.745x10® 

0.05 

♦ 

ay 

13.5x10® 

6.75x10^ 

0.05 

_♦ 

H 

60.0x10® 

3.0x10® 

0.05 

♦ 

it 

0‘ 

O 

o 

— 

* 

6 

7.95x10“^ 

— 

0.05 


♦(Ty indicates yield stress, H indicates hardening modulus, 0 indicates 

ply orientation angle, and 6 indicates ply thickness. 

The values of the a- ■ constants in the yield criterion are: 

|a^^ = 0.001 , §a 22 = 10 , |aj 2 = -0-01 


5^44 




55 




66 


= 1.9 


The hardening model used was Y(a) — H[a + 


jfYjA]l/A 

H 


A = 5.8 


Units are in psi and inches where appropriate. 



34 


UWIMTE 



t 

••SYMTHESIS' 


f UMIHATE 
• THEORY 



PLY 


\ 


MATERIAL PROPERTIES 
P W, T, t) 



MONTE CAmO 

\ SAMPLING 

/ 

f 


LAMINATE 

THEORY 


i 

“DECOMPOSmON" 



PLY 


/ 


X 

COMPOSITE 

MtCROMECHANICS 

THEORY 



CONSTITUENTS 




COMPOSITE 

MICROMECHANICS 

THEORY 


Figure 1. Probabilistic integrated multi-scale metal 
matrix composite analysis. 
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Figure 2. Multi -factor constituent material model and 
micromechanics subcell. 
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fa) Case 1 histogram (narrow range) 
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(b) Case 1 cumulative distribution 
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(c) Case 2 histogram (wide range) 
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(d) Case 2 cumulative distribution 


Figure 3. Histograms and cumulative distribution curves for 
in-plane shear modulus 6,^2 
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Gr-Cu CONFIDENCE INTERVALS 



Figure A. Longitudinal tensile strength 
with perturbed shape parameter of fiber 
strength. 


Gr-Cu CONFIDENCE INTERVALS 



FIBER VOLUME RATIO 


Figure 5. Longitudinal compressive strength 
with perturbed shape parameter of matrix 
strength. 
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Gr-Cu CONFIDENCE INTERVALS 



Figure 6. Transverse tensile strength with 
perturbed shape parameter of matrix strength. 


Figure 7. Transverse compressive strength 
with perturbed shape parameter of matrix 
strength. 
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Figure 8. Longitudinal tensile strength 
with perturbed COV of fiber angle. 
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Figure 10. Transverse tensile strength 
with perturbed COV of fiber angle. 
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Figure 9. Longitudinal compressive strength 
with perturbed COV of fiber angle. 
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Figure H* Transverse compressive strength 
with perturbed COV of fiber angle. 
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Figure 12. Longitudinal tensile strength 
with perturbed shape parameter of fiber 
strength; power law inactive. 
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Figure 14. Longitudinal tensile strength 
with perturbed COV of fiber melting temperature; 
temperature power term active. 
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Figure 13. Longitudinal tensile strength 
with perturbed COV of fiber strength temperature 
exponent n; temperature power term active. 
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Figure 15. Longitudinal tensile strength 
with perturbed shape parameter of fiber 
strength; stress power term active. 
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Figure 16. Longitudinal tensile strength 
with perturbed COV of matrix modulus; power 
law inactive. 
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Figure 18. Longitudinal tensile strength 

with perturbed COV of matrix melting temperature; 

temperature power term active. 
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Figure 17. Longitudinal tensile strength 
with perturbed COV of matrix modulus temperature 
exponent n; temperature power term active. 


SCS-6 Til 5 CONFIDENCE INTERVALS 



Figure 19. Longitudinal tensile strength 
with perturbed shape parameter of mechanical 
strength; mechanical cycle power term active. 
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Figure 20. ARALL laminate layup anc geometry. 
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Figure 21. ARALL tesnsion test experimental and 
analytical comparisons. 
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Figure 22. Finite element mesh, loading and boundary 
conaitions of a tension specimen with 

a hole. 
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Figure 23. Mean and standard deviation of normal strain 

at the hole edge {point A) versus load for the 
case where the fibers are at 90“ to the load. 
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Figure 24. Mean and standard deviation of normal strain at 
the hole edge (point A) versus load for the case 
v/here fibers are aligned with the load. 



Figure 25. Analytical and experimental comparison of the 
longitudinal strain along a line 90° to the 
loading for the Boron/Aluminum laminate. 
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Figure 26.Di stn buti cn of the longitudinal strain along 
a ;:re 90' to the loading for the Boron/Al 
laminate tension specimen with hole. 
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Figure 27. Mean and standard deviation of normal strain 

at tne hole edge (point A) versus load for the 
Boron/Al umi num tension specimen with hole. 
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